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Abstract — It is described, explicitly, how a popular, 
commercially-available software package for solving partial- 
differential-equations (PDEs), as based on the finite-element 
method (FEM), can be configured to calculate the frequencies 
and fields of the whispering-gallery (WG) modes of axisymmetric 
dielectric resonators. The approach is traceable; it exploits 
the PDE-solver's ability to accept the definition of solutions 
to Maxwell's equations in so-called 'weak form'. Associated 
expressions and methods for estimating a WG mode's volume, 
filling factor(s) and, in the case of closed(open) resonators, its 
wall(radiation) loss, are provided. As no transverse approxi- 
mation is imposed, the approach remains accurate even for so- 
called quasi-TM and -TE modes of low, finite azimuthal mode 
order. The approach's generality and utility are demonstrated by 
modeling several non-trivial structures: (i) two different optical 
microcavities [one toroidal made of silica, the other an AlGaAs 
microdisk]; (ii) a 3rd-order microwave Bragg cavity containing 
alumina layers; (iii) two different cryogenic sapphire X-band 
microwave resonators. By fitting one of the latter to a set 
of measured resonance frequencies, the dielectric constants of 
sapphire at liquid-helium temperature have been estimated. 



I. Introduction 

EXPERIMENTAL data are related to physical laws, ex- 
pressed as equations, through models. To determine the 
either fundamental, phenomenological, or 'materials' con- 
stants that the model's equations include, the model must first 
be solved, and explicitly so, to allow the fitting of its constants, 
through (Bayesian) regression, to the experimental data. The 
inaccurate solution of a model can sometimes contribute 
significantly to, if not wholly dominate, the fitted values' 
uncertainties. Improvements in the accuracies of solutions can 
alone motivate the (re-)determinations of constants from extant 
(i.e. 'old') experimental data. Indeed, the method of solution 
presented in this paper's section II is subsequently exploited, 
in section VI, to determine the values of certain dielectric 
constants from the frequencies at which an electromagnetic 
resonator was found to resonate experimentally; here, the 
'modeling errors' dominate over other uncertainties. 

Once all relevant physical constants are known to sufficient 
accuracy, a model's solution can be exploited in the reverse 
sense to simulate as-yet unrealized experimental embodiments. 
Simulation enables the properties of a proposed embodiment 
to be predicted and, thus, through modifications, for its per- 
formance (with respect to a given application) to be optimized 
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without, or at least at reduced, experimental effort. Suffi- 
cient accuracy in the model's solution itself is, again, vital. 
Though analytical models can adequately treat certain highly 
symmetrical structures, the sufficiently accurate solution of 
less symmetrical (though more practical) structures typically 
requires automated numerical computation -as implemented 
on a digital computer. Here, the model's structure must first 
be represented in some 'electronic' format. Then, the physical 
equations, often including sets of partial-differential ones, 
are encoded and solved for the boundary conditions and 
constitutive relations that the structure implies. Though these 
two tasks can be implemented by hand-coding in low-level 
computer languages, highly developed commercial software 
packages now exists to facilitate both: (i) computer-aided- 
design (CAD) tools and (ii) partial-differential-equation (PDE) 
solvers, respectively. Many, though by no means all [1], [2], 
[3], of the latter are based on the finite-element method (FEM) 
[4], which can readily accept CAD-defined structures. 

Furthermore, various packages now integrate (i) and (ii) 
into complete computer-modeling environments, e.g. 'ECAD' 
for simulating electromagnetic systems. These environments 
sport various additional features for accelerating the defini- 
tion of models and for facilitating the display and analysis 
('post-processing') of solutions; they also impose standardized 
formats and provide ('house-keeping') tools to assist in the 
maintenance, sharing and documentation of a model -so that 
others can subsequently benefit from, and build upon, the 
original model-developer's effort. Compared to the laborious 
coding up and piecing together of the equivalent software by 
hand (e.g. as straight MATLAB or Fortran code, making calls 
to optimized 'canned' matrix eigensolvers), the use of such 
environments, despite their costs and limitations, is attractive. 

A problem associated with the inclusion of a complex model 
into the determination of a constant, where the model is solved 
via a piece of commercial 'black-box' software, or through 
a proprietary 'in-house' service, is that the determination 
may thus cease to be traceable. Significant effects (or 'un- 
documented features') imparted by the modeling/simulation 
process may become difficult if not impossible to isolate, 
understand, or quantify. With regard to traceability, both the 
model's definition and its chosen method of solution must re- 
main amenable to explicit representation, thus communication, 
thus external scrutiny. Convenience and/or efficiency demand, 
furthermore, that this representation be as concise and elegant 
as possible -with no ambiguities. 
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A. Whispering-gallery-mode resonators 

Electromagnetic structures that support whispering-gallery 
(WG) modes are technologically important because of the 
advantageous properties that these modes exhibit in terms 
of spatial compactness, frequency control (either stability or 
agility) and mode quality (Q) factor. Explicit examples will 
be presented and/or referenced in due course. Compared to 
the abrupt retro-reflection of an electromagnetic wave at the 
surface of a mirror, the continuous bending of the same by a 
whispering-gallery waveguide is an alternative, and at that a 
still relatively underdeveloped one, which is opening up new 
applications. Here, one is often interested in a 'resonator' 
where the wave's trajectory closes back on (and the wave 
thereby interferes with) itself. Though elliptical [5], helical, 
or even more complex bending trajectories [6] can be (and 
have been) envisaged in association with the various mor- 
phologies of electromagnetic waveguide/resonator that support 
WG waves, the author restricts himself here to the study of 
the simplest, and to-date most popular, class of WG-mode 
resonator: that where the electromagnetic wave's trajectory 
is a plane circle (thus constant radius of curvature) and the 
electromagnetic structure supporting it is axisymmetric (and 
coaxial with respect to the said circle/WG mode). 1 Within 
this class, a convex dielectric:vacuum boundary is often the 
curved interface of choice for guiding/confining the whispering 
gallery mode around in an circle. The method presented 
below can, however, also be employed to simulate WG-modes 
that are guided by a concave dielectric:metal boundaries. In 
general then, one considers an axisymmetric toroidal volume, 
whose cross-section in a (it matters not which) radial-axial 
plane comprises regions of dielectric (voids correspond to the 
dielectric vacuum) that are bounded (either externally or from 
within) by metal surfaces; see FIG. 1. 

Despite the breadth and technological allure of this class of 
WG-mode resonator, it is the author's understanding that most 
if not all commercial (ECAD) packages available at the time of 
writing (early to mid. 2006) suffer from a rather unfortunate 
'blind spot' when it comes to calculating, efficiently (hence 
accurately), the whispering-gallery modes (with plane circular 
trajectories) that such axisymmetric resonator's support. The 
popular MAFIA/CST package [7], with which the author is 
familiar, is a case in point: As has also been experienced 
by Basu et al [8], and no doubt others, one simply cannot 
configure the software to take advantage of the WG modes' 
fl/?n'on'-known azimuthal dependence, viz. exp(iM0), where 
M is a positive integer known as the mode's azimuthal mode 
order, and <fi is the azimuthal coordinate. Though frequencies 
and field-patterns can be obtained (at least for WG modes 
of low azimuthal mode order), the computationally advanta- 
geous reduction of the problem from 3D to a 2D that the 
rotational symmetries of the resonator and its solutions allow 
is, consequentially, precluded; and the ability to simulate high- 
order whispering-gallery modes with sufficient accuracy for 
metrological purposes is, exasperatingly, lost. About the best 

'it is acknowledged that even axisymmetric (3D) resonators can support 
'spooling' -helical WG modes that do not lie in a plane [6]; the analysis of 
such exotica lies outside the scope of this paper. 



one can do is to simulate a 'wedge' [over an azimuthal domain 
A<fi — n/(2M) wide] between radial electric and magnetic 
walls [9]. 

From ECAD to 'omni-CAD': Adding titillation to the exas- 
peration, several commercial packages [10], [11] based on 
the FEM method are now beginning to offer true 'multi- 
physics' capabilities: not only can one separately model a 
resonator's electromagnetic response, its mechanical response, 
its thermal response, all based on a common, defined- 
once-and-for-all geometric structure, one can furthermore cou- 
ple/ 'extrude '/integrate these responses to model non-linear 
and/or parametric effects. These effects include (as illustrative 
examples): (i) the electromagnetic heating of a resonator's 
lossy dielectrics and/or it resistively conducting inner surfaces 
(thus shifts in the frequencies of the resonator's electromag- 
netic modes), and -even- (ii) 'mechanical-Kerr' instabili- 
ties/oscillations associated with the mechanical deformation of 
the resonator's components due to radiation pressure [12], as 
exerted by a driven electromagnetic mode. This nirvana of pre- 
dictive (+ deductive) capability is, for WG-mode resonators, 
in view of the alluring applications associated with their 
nonlinear and/or parametric effects, a particularly tantalizing 
destination -if only one could appropriately configure the (in- 
the-first-place sufficiently configurable) software to get there. 
This paper provides a single, though -one might claim- a quite 
fundamental, generic, and enabling, step on the long march 
there to. 

B. Brief, selected history of WG-mode simulation 

The analysis and modeling of the whispering-gallery modes 
of electromagnetic resonators, at both optical and microwave 
frequencies, continues to support and guide experimen- 
tal endeavor [13], [14], [15], [16]. A brief and far-from- 
comprehensive survey of the different methods used to im- 
plement these simulations to date, with a strong selective 
bias towards those that have been applied to the study of 
microwave dielectric-ring resonators, is provided here. The 
author encourages the reader to consult the earlier works that 
are referenced within the papers cited below. 

Based on 'separating the variables' (SV), textbooks [17], 
[18] provided expressions for the whispering-gallery modes 
that are supported by whole dielectric parallelopipeds, right- 
cylinders, and spheres, or dielectric layers and shells exhibiting 
the same symmetries thereof, where the dielectric volume is 
enclosed by electric and/or magnetic walls. Illustrating the 
genre is Wilson et afs handy study of the transverse electric 
(TE) and transverse-magnetic (TM) modes of right-cylindrical 
metal cavities [19], which was in fact used by the author to 
validate the weak-form expressions described in sub-section 
II-A below in the early stages of his work. 

By 'mode-matching' (MM) these SV-solutions across 
boundaries, WG-mode solutions for composite axisymmetric 
structures, such as a dielectric right cylinder, surrounded 
immediately by a void, enclosed within a (coaxial) right- 
cylindrical metallic jacket {i.e., a so-called 'can'), can be 
derived [20], [21]. These solutions, with their associated dis- 
crete/integer indices (related to symmetries), provide a nomen- 
clature [22] for classifying WG modes. This nomenclature can 
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be re-used to sort and label the lower-order WG modes of 
less symmetrical though structurally similar resonators, where 
these modes can only be calculated 'blindly' through other, 
more numerical methods. Mode-matching by taking linear 
combinations of several/many -as opposed to just a few- basis 
functions can increase the 'fit' hence accuracy of the MM-SV 
method and/or allow it to be extended it to the treatment of 
deformed structures. 

In view of the limited computational effort that these semi- 
analytical SV-based methods demand, remarkable accuracies 
can be achieved, especially when the most 'sympathetic' basis 
functions are deployed. For many resonators of interest, how- 
ever, and the composite axisymmetric structures mentioned in 
the previous paragraph are a case in point, it is not possible, 
in the MM-SV method (as based on a finite set of basis 
functions), to simultaneously match all components of the 
electric and magnetic fields across all boundaries [15] -to 
do so would, after all, amount to an exact solution! With a 
small, finite basis set, the 'transverse' (or 'axial-polarization') 
approximation, that tolerates a mis-match of 'minor' field 
components, whilst consistently matching the major ones, 
is uncontrolled. Though extensions to mode-matching that 
capture spatially non-uniform polarization can be constructed 
[3], the MM-SV method in general needs to be validated, for a 
given shape of resonator and mode, through comparison with 
(more exact) solutions supplied by other methods. 

For a complete, accurate solution of Maxwell's equations, 
one must generally resort to wholly numerical methods, of 
which there are several relevant classes and variants. Apart 
from the finite-element method (FEM) itself (considered in 
more detail further on), the most developed and thus most 
immediately exploitable alternatives include (given here for 
reference -not considered in any greater detail): (i) the Ritz- 
Rayleigh variational or 'moment' methods [23], [24], [25], 
(ii) the finite difference time domain method (FDTD) [1], 
[26], and (iii) the boundary-integral [2] or bounday-element 
methods (BEM, including FEM-based hybridizations thereof 
[27]). Zienkiewicz and Taylor [4], though nominally dedicated 
to FEM, provide a taxonomy {viz. table 3.2 loc. cit.) covering 
most of these methods, which reveals certain commonalities 
between them. It is remarked here, for example, that FDTD 
may be regarded as a variant of the FEM employing local, 
discontinuous shape functions. 2 In conjunction with these 
generalities, it is worth re-iterating here that the core for- 
mulation presented in this paper (viz. equations 8 through 
23) can exploit any PDE-solver (e.g. a moment-method-based 
one) capable of accepting/intepreting weak-form statements. 
Though a FEM-based solver (viz. COMSOL/FEMLAB) was 
indeed used to provide the examples presented in section V 
below, this article's formulation is not, per se, wedded to FEM. 

Though the finite-element method can solve for all field 
components (both major and minor), the explicit, direct state- 
ment of the required set of a coupled partial differential equa- 

2 It is also worth remarking that, for resonators comprising just a few, large 
domains of uniform dielectric, then the boundary-integral methods (based on 
Green functions), which -in a nutshell- exploit such uniformity to reduce 
the problem's dimensionality by one, will generally be more computationally 
efficient than FEM. 



tions (i.e. Maxwell's equations for WG-mode electromagnetic 
resonators) in component form, suitable for the insertion into 
a standard commercial FEM/ECAD software package, can 
be extremely onerous -if not absolutely ruled out by the 
software's lack of configurability. This is why the majority 
of these packages already include pre-defined 'applications 
modes', 'macros' or 'wizards' for solving electromagnetic 
problems (for particular geometries). To simplify the problem 
to that of a single (scalar) PDE, one can again [cf. the 
SV-MM method(s) already discussed above] invoke the so- 
called transverse (axial-polarization) approximation, where the 
resonator's either magnetic or electric field is forced to lie 
everywhere parallel to the resonator's axis of rotational sym- 
metry; figure B.l of reference [28] displays this approximation 
most pedagogically. Investigations based on this 'transverse- 
FEM' approach have been reported in several recent works [8], 
[28]. 3 Though it can provide indicative trends and quantitative 
results, which might well be sufficiently accurate and/or robust 
for the calculation's intended purpose (in view of even less 
well controlled experimental parameters), the uncontrolled 
nature of the approximation that transverse-FEM incorporates 
is again far from ideal. The careful physicist, or metrologist, 
is (again) compelled to justify its validity, for a given res- 
onator and mode, through comparison with either non-trivial 
analytical solutions, where they exist, or 'brute force' (3D) 
numerical computation [8]. It is this paper's principal claim 
that, through only a modicum of extra effort, the transverse 
approximation, and its associated onerous validations [or (else) 
lingering doubts], can be wholly obviated. 

The application of the finite-element method to the solving 
of Maxwell's equations has a history [29], and is now very 
much an industry [7], [10], [11]. Zienkiewicz and Taylor 
[4] supply FEM's theoretical underpinnings, in particular an 
erudite account of Galerkin's method of so-called 'weighted 
residuals'. A pervasive, and often quite debilitating problem 
that besets the direct/'naive' applications of FEM to the PDEs 
that are Maxwell's equations is the generation of (a great 
many) spurious solutions [30], [31], associated with the local 
gauge invariance, or 'null space' [31], which is a (hidden) 
feature of these PDEs (in particular its 'curl' operators). 
At least two research groups have nevertheless successfully 
developed software tools for calculating the WG modes of 
axisymmetric dielectric resonators, where these tools (i) solve 
for all field components (i.e. no transverse approximation is 
invoked), (ii) are 2D (and thus numerically efficient) and (iii) 
effectively suppress spurious solutions (without any insidi- 
ous, detrimental side-effects) [30], [23], [24], [32], [33]. The 
method reported below sports these same three attributes. With 
regard to (iii), the approach adopted by Auborg et al was 
to use different finite elements (viz. a mixture of 'Nedelec' 
and 'Lagrange' -both 2nd order) for different components 
of the electric and magnetic fields; Osegueda et al [32], on 
the other hand, used a so-called 'penalty term' to suppress 
(spurious) divergence of the magnetic field. Stripping away 
all of its motivating remarks, applications and illustrations, the 

3 Srinivasan et al [14] in contrast employ a 'full-vectorial model' though 
they do not explicitly define what this is. 
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principal function of this paper is to convert (one might say 
'extract') the method encoded by Osegueda et al's 'CYRES' 
2D FEM software package [32] into explicit 'weak-form' 
expressions, that can be directly and openly ported to any PDE 
solver (most notably COMSOL/FEMLAB) capable of accept- 
ing such. These weak-form expressions are wholly equivalent 
to the Maxwellian PDEs from which they are derived. But, 
being scalar (tensorially-contracted), they are considerably less 
onerous to represent and communicate than the vectorial PDEs 
themselves. The author hopes that, by stating/popularizing the 
problem so explicitly in this paper through the lingua franca of 
weak form, the means to model, both accurately and traceably, 
the whispering-gallery modes of axisymmetric resonators will 
thereby be made accessible to any competent engineer or 
physicist in need of such a means -'off the shelf, as opposed 
to it remaining a strictly '/n-house' (and thus rather less open 
and traceable) capability retained by specialists. 

II. Method of solution 

A. Weak forms 

Scope: The types of ideal resonator that fall under the scope 
of the analysis presented immediately below are those that 
comprise volumes of lossless dielectric space bounded by 
a combination of perfect (thus also loss-less) electric or 
magnetic walls -see Fig. 1 (though note that the restriction 
to resonators of axisymmetric form needs only to be invoked 
at the start of subsection II-B). As discussed in subsections Ill- 




Fig. 1. Generic axisymmetric resonator in cross-section (medial half- 
plane). A dielectric region enclosed by an electric wall El is subdivided 
into subdomains Dl, D2 and D3, none, one or more of which could be 
free space. D2 and D3 are bounded internally by electric walls E2, E2' 
and E3. The resonator's mirror symmetry through the horizontal (equatorial) 
plane, as indicated by the dashed line Ml, allows an additional either electric 
or magnetic wall to be placed on it and, thereupon, only one half of the 
resonator's domain (either its upper or lower half) need be simulated. 

D through IV-B below, a dissipative (open) resonator's finite 
fractional energy loss per cycle, hence its sub-infinite Q-factor, 
can be estimated, and with often perfectly sufficient accuracy 
with regard to applications, from the solution of an equivalent 
loss-less (closed) resonator. The resonator's dielectric space 
will, in general, comprise both voids (i.e. the free space of a 
vacuum) and space filled by sufficiently 'good' (i.e. low-loss) 



dielectric material(s). A model's electric walls will translate, 
in embodiments, to metallic surfaces whose conductivities are 
sufficiently good to be treated as such (section III-D quantifies 
the loss caused by the metallic wall of a particular resonator's 
can). The (relative) permeability of real magnetic materials is 
seldom high enough for a wall made from such to be regarded, 
without deleterious approximation, as a (perfect) magnetic 
one. When modeling resonators whose forms exhibit reflection 
symmetries, such that the magnetic and electric fields of their 
solutions exhibit either symmetry or antisymmetry through 
mirror planes, perfect magnetic and electric walls can be 
advantageously imposed on the model's mirror plane to solve 
for particular 'sectors' of solutions. 

The electromagnetic field within the dielectric volumes of 
the resonator obeys Maxwell's equations [34], [18], as they 
are applied to continuous, macroscopic media [35]. Thus, 
on the assumption that the resonator's constituent dielectric 
elements have negligible (or at least the same) magnetic 
susceptibility (hence permeability), the magnetic field strength 
H is continuous across all interfaces between them. 4 This 
property makes it advantageous to solve for H (or the magnetic 
flux density B = /iH related to it by a constant global magnetic 
permeability f£), as opposed to the electric field strength E (or 
displacement D). Upon substituting one of Maxwell's 'curl' 
equations equations into the another, the problem reduces to 
that of solving the (modified) vector Helmholtz equation 

V x (e^ 1 V x H) - aV(V • H) + C - 2 d 2 tL/dt 2 = (1) 

subject to appropriate boundary conditions (read on), where 
c is the speed of light. Here, e _1 denotes the inverse relative 
permittivity tensor; one assumes that the resonator's dielectric 
elements are linear, such that it is a (tensorial) constant - 
i.e. independent of field strength. The middle (penalty) term 
on the left-hand side of equation 1 is the same as that used 
by Osegueda et al [32]; it functions to suppress/reveal spu- 
rious modes in the finite-element simulation 5 ; the constant a 
controls this term's weight with respect to its Maxwellian neig- 
bours. The penalty term's insertion into the above Helmholtz 
equation is wholly permissible since, for every true solution 
of Maxwell's equation, it must exactly vanish (everywhere): 
the magnetic flux density B, hence (for non- or isotropically- 
magnetic media) H = B/fi, is required to be free of divergence 
-assuming no magnetic monopoles lurk inside the resonator. 

Reference [18] (particular section 1.3 thereof) supplies a 
primer on the electromagnetic boundary conditions discussed 
forthwith. Assuming that the resonator's bounding electric 
walls are perfect in the sense of having (effectively) infinite 
surface conductivity, the magnetic flux density at any point 
on each such wall is required to satisfy B • n = 0, where 
n denotes the wall's surface normal vector. Providing the 

4 The method described in this paper could be extended (still within the 
configurability of the PDE-solver) to the analysis of resonators containing 
different magnetic dielectrics (such as YIG, and as would exhibit differing 
permeabilities) by way of 'coupling variables' -to transform H (or B) across 
internal boundaries between regions of different permeability. 

5 It is mentioned here that the author briefly experimented in COMSOL 
with mixed ('Nedelec' plus 'Lagrange') finite elements [23] but found that the 
above penalty term, in conjunction with 2nd-order Lagrange finite elements 
applied to all three components of H, gave wholly satisfactory results. 
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magnetic permeability/susceptibility of the dielectric medium 
bounded by the electric wall is not anisotropic, this condition 
is equivalent to 

H • n = 0. (2) 

The electric field strength at the electric wall is required to 
obey 

E x n = 0; (3) 

these two equations ensure that the magnetic (electric) field 
strength is oriented tangential (normal) to the electric wall. 
As is pointed out in reference [32], equation 3 is a so-called 
'natural' (or, synonymously, a 'naturally satisfied') boundary 
condition within the finite-element method -see ref. [4]. 

When the resonator's form (hence those of its solutions) 
exhibits one or more symmetries, it is often advantageous 
(for reasons of computational efficiency) to solve only for 
a symmetry-reduced portion or 'sector' of the full resonator, 
where this sector is bounded by either (real or virtual) electric 
walls or (virtual) magnetic walls, or both. The boundary 
conditions corresponding to a perfect magnetic wall (dual to 
the those for an electric wall) are 

D • n = 0, (4) 

and 

H x n = 0; (5) 

these two equations ensure that the electric displacement 
(magnetic field strength) is oriented tangential (normal) to the 
magnetic wall. Again, the latter equation is naturally satisfied. 

One now invokes Galerkin's method of weighted residuals; 
reference [4] explains the fundamentals here; reference [31] 
provides an analogous treatment when solving for the electric 
field strength (E). Both sides of equation 1 are multiplied 
(scalar-product contraction) by the complex conjugate of a 
'test' magnetic field strength H , then integrated over the 
dielectric resonator's interior volumes. Upon expanding the 
permittivity-modified 'curl of a curl' operator (to extract a 
similarly modified Laplacian operator), then integrating by 
parts (spatially), then disposing of surface terms through the 
electric- or magnetic-wall boundary conditions stated above, 
one arrives (equivalent to equation (2) of reference [32]) at 

/ [(V x H*) - (V x H) — 

Jv e 

a(V ■ H*)(V • H) + c~ 2 H* • d 2 tL/dt 2 } dV = 0, (6) 

where ' f y ' denotes the volume integral over the resonator's 
interior space (or sector thereof) and '-/e' denotes a con- 
traction weighted by inverse relative permittivities. The three 
terms appearing in the integrand correspond directly to the 
three weak-form terms required to define an appropriate finite- 
element model within the PDE solver. 

Assuming that the physical dimensions and electromagnetic 
properties of the resonator's components are temporally invari- 
ant (or at least 'quasi-static'), one seeks harmonic or 'modal' 
solutions: H(r;t) = H(r)exp(— i2wft), where r is the vector 
of spatial position, t the time, and / the mode's resonance 
frequency. The last, 'temporal' term in equation 6's integrand 



can thereupon be re-expressed as — (c/) 2 H(r)* • H(r), where 
c = 2-k/c and c is the speed of light. This re-expression(and, 
with respect to Spillane et al's work, using exactly the same 
FEM software platform) reveals the integrand's complete dual 
symmetry between H and H. 

B. Axisymmetric resonators 

One now restricts the scope of the analysis to resonators 
whose interiors and bounding surfaces are electromagnetically 
axisymmetric (henceforth referred to simply as 'axisymmet- 
ric resonators') where a system of cylindrical coordinates 
is aligned with respect to the resonator's axis of rotational 
symmetry. This system's three components are {x, <j>, y} = 
{'rad(ial)', 'azi(muthal)', 'axi(ial)'} 6 . One wishes to calcu- 
late the resonance frequencies and field patterns of the res- 
onator's (standard) whispering-gallery (WG) modes, whose 
phase varies as exp(iM0), where M — {0, 1, 2, ...} is the WG 
mode's azimuthal mode order. Note that the method does not 
require M to be large (i.e., it is not an 'asymptotic' method); 
even modal solutions that are themselves axisymmetric, cor- 
responding to M = 0, such as the one shown in Fig. 6(b), 
can be calculated. Viewed as a three-component vector field 
over a (for the moment) three-dimensional space, the time- 
independent part of the magnetic field strength now takes the 
form 

H(r) = e lM * { H x {x, y),i H+{x, y),H«{x, y) } (7) 

where an T (= 1)) has been inserted into the field's 
azimuthal component so as to allow, in subsequent solutions, 
all three component amplitudes {H x ,H' t ',H y } to each be 
expressible as a real amplitude multiplied by a common 
complex phase factor. The relative permittivity tensor of 
an axisymmetric dielectric material is diagonal with entries 
(running down the diagonal) ed; ag . = {e±,e±,e^}, where ey 
is the material's relative permittivity in the axial direction and 
e± its relative permittivity in the plane spanned by it radial 
and azimuthal directions. 

It now remains only to substitute equation 7 into equation 
6's integrand and express the three terms composing the 
latter's integrand in terms of the magnetic field strength's com- 
ponents (and their spatial/temporal derivatives); textbooks pro- 
vide the required explicit expressions for the vector differential 
operators in cylindrical coordinates [36], [17], [18]. A radial 
factor, x, is included here from the volume integral's measure: 
dV = 2n x dx d(f> dy (the factor of 2n here is uniformly, thus 
inconsequentially, dropped from all three expressions below.) 
These requisite expansions are presented here in compact 
mathematically notation; their line-text (i.e. with no super- or 
sub-scripts, hence considerably more verbose) equivalents, in 
forms suitable for direct cut-and-paste injection into a popu- 
lar PDE-solver (viz. COMSOL/FEMLAB) are available as a 
separate Appendix' to this paper [37]. The first, 'Laplacian' 
weak term is given by 

(V xH)-(VxH) = (^|B + xC)/(e ± e\\), (8) 

6 'x' and 'j/', instead of the more conventional V and 'z', are (regrettably) 
used to represent radial and axial coordinates/components, respectively, so as 
to comply with COMSOL/FEMLAB's standard (2D) naming conventions. 
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where 

A = { e^H'f'H 4 ' - M{B*B X + B*B X ) + M 2 B X B X )\ 

+ e ll M 2 H y H y }, (9) 

B = e±[H^(H^ - MH X ) + H^(H^ - MB X )] 

-e\\M{H y H* + H y H%), (10) 

C={e^HtHi 

+ e„ \{H y - B x )(B y - B x ) + H$H$] }, (11) 

where H% denotes the partial derivative of H * (the aziumthal 
component of the magnetic field strength) with respect to 
x (the radial component of displacement), etc.. Here, the 
individual factors and terms have been ordered and grouped so 
as to display the dual symmetry. Similarly, the weak penalty 
term is given by 

a(V • H*)(V • H) = a{— + E + x F}, (12) 

x 

where 

D = H X H X - M{B*B X + B*B X ) + M 2 H 4 'H' t ', (13) 



E=(B X + B y )(B x - MB*) 



+ (B X - MB*)(B X +B y ), 



F^(B x + B y )(B x +B y ). 



(14) 
(15) 



And, finally, the temporal weak-form ('dweak') term is given 

by 

H-d 2 U/d 2 t = c- 2 x (H x H? t +H' t 'H? t +H y H? t ) 

= -c 2 f 2 x (B X B X + B^B* + B V B V ), (16) 

where B x t denotes the double partial derivative of B x 
w.r.t. time, etc.. What is crucial is that none of the terms on 
the right-hand sides of equations 8 through 16 depend on the 
azimuthal coordinate <f>; the problem has been reduced from 
3D to 2D. 

C. Axisymmetric boundary conditions 

An axisymmetric boundary's unit normal in cylindrical 
components can be expressed as {n x ,0,n y } -note vanishing 
azimuthal component. The full electric-wall boundary condi- 
tions, in cylindrical components, are as follows: H n = 
gives 

H x n x + H y n v = 0, (17) 



and E x n = includes both 



B x -B y = 



(18) 



and 



When the dielectric permittivity of the medium bounded by 
the electric wall is isotropic (which is often the case in 
embodiments), the last condition reduces to 

(H* - B X M + H$x)n x - (B y M - B*x)n y = 0. (20) 

The full magnetic-wall boundary conditions, in cylindrical 
components, are as follows: the condition D n = gives 

(B y M - H+x)n x + (B* - B X M + Bfx)n y = 0, (21) 

and the condition H x n = includes both 



and 



B y n x - B x n y = 



H* = 0. 



(22) 



(23) 



One observes that the transformation {n x — > — n y , n y — ► n x }, 
connects equations 17 and 22, and equations 20 with 21. Ex- 
plicit PDE-solver-ready equivalents of 17 through 23 are stated 
in this paper's auxiliary Appendix [37]. The above weak-form 
expressions and boundary conditions, viz. equations 8 through 
23 are the key enabling results of this paper: once inserted 
into a PDE-solver, the WG modes of axisymmetric dielectric 
resonators can readily be calculated, as is demonstrated for 
particular embodiments in section V below. 

III. Postprocessing of solutions 

Having determined, for each mode, its frequency and all 
three components of its magnetic field strength H as functions 
of position, other relevant fields and parameters can be derived 
from them. 

A. Remaining Maxwellian fields 

Straightaway, the magnetic flux density B = H/fi; here, 
as stated in subsection II-A above -but see also footnote 
4, the magnetic permeability \i is assumed to be a scalar 
constant, independent of position. [And for each of the res- 
onators considered in section V, n = [1q everywhere -to 
an adequate approximation.] As no real ('non-displacement') 
current flows within a dielectric, V x H(t) = dD(i)/dt, thus 
D = -i(27r/) _1 V x K(t). And, E = e -1 D, where e" 1 is the 
(diagonal) inverse permittivity tensor, as already discussed in 
connection with equation 6 above. 

B. Mode volume and filling factor 

Accepting various caveats (most fundamentally, the problem 
of mode-volume divergence -see footnote 7; and inconsistent 
definitions between different authors ...) as addressed by Kip- 
penberg [28], the volume of a mode in a dielectric resonator 
is here defined as [14] 



^modc 



Jjj h _ s e\E\ 2 dW 
max[e|E| 2 ] 



(24) 



e ± (H' l '-H x M + H$x)n x -e ll (H y M-H*x)ny = 0. (19) 



where max [...], denotes the maximum value of its functional 
argument, and J J J h _ s ...dV denotes integration over and 
around the mode's 'bright spot' -where its electromagnetic 
field energy is concentrated. 
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C. Filling factor 

The resonator's electric filling factor, for a given mode, 
a given dielectric piece/material, dicl., and a given field 
direction, (dir. € {radial, azimuthal, axial}), is defined as 

J7/e|E|»dV ' ^ 

where J J J diol ...dV denotes integration only over those do- 
mains composed of the dielectric in question and pol. = {_L, ||} 
for dir. = {radial or azimuthal, axial}. The numerators and 
denominators of equations 24 and 25 can be readily evaluated 
using the PDE-solver's post-processing features. 

D. Finite Qs and wall losses 

So far, the model resonator as per Fig. 1 has been treated 
as a wholly loss-less one: its modal solutions have infinite 
Qs or, equivalently, the (otherwise complex) frequencies of 
these solutions are purely real or, equivalently, the solutions' 
oscillatory electromagnetic fields persist indefinitely. No en- 
ergy is dissipated by the dielectric material(s) included within 
the resonator (their dielectric loss tangents are presumed to 
be zero); none is lost through radiation -since the resonator's 
bounding perfect electric/magnetic walls allow none to escape; 
and, being perfect, the current induced within the wall causes 
no resistive dissipation. 

Real resonators, on the other hand, are subject to one or 
several dissipative processes, i.e. 'losses', that render the Qs 
of resonances finite. This subsection provides an expression 
for the rate of a resonator's 'wall loss'; section IV goes 
on to provide bounds on the rate of an open resonator's 
'radiation loss'. Such estimates are important since Q values 
are directly measurable in experiments and, furthermore, often 
determine viability and/or performance in applications. The 
approach taken here is to build upon (via perturbation theory 
and/or 'induction') the loss-less model, as it has already 
been formulated in section II -as opposed to constructing the 
additional machinery required to model resonators with lossy 
materials either placed within or clad about them. 

Preliminaries: The space that a WG-mode occupies can be 
broadly divided into three regions: (i) the mode's 'near-field', 
which includes its bright spot(s) -where the modes e.m. energy 
density is greatest, (ii) an 'evanescent zone', lying around 
the near-field, where the mode's field amplitudes (and energy 
density) decay exponentially with distance r away from its 
bright spot and (iii) a (notionally infinite) 'radiation zone', 
lying outside of the evanescent zone (a 'cusp' can separate the 
two), where the field-amplitudes decay as ~ 1/r. If a compact, 
closed resonator of high-Q is sought, its electric (i.e. metal) 
wall should be placed sufficiently far from the WG-mode's 
bright spot in the exponentially decaying zone (ii), but, for 
reasons of compactness, no further out than is necessary. 

A brief word of warning: As the Q of an experimental WG- 
resonator can be exceedingly high (> 10 9 ), the amplitude of 
the electromagnetic field where dominant losses occur (their 
rates will depend on the amplitude) can be many orders of 
magnitude lower than the field's maximum (or maxima) at 
the WG-mode's bright spot(s). The hardware and software 



employed to generate WG-modal solutions must thus be able 
to cope with such a dynamic range, lest significant numerical 
errors creep into the predicted (loss-rate-determining) ampli- 
tudes where the WG-resonator's supported mode is faint. In 
practice, this means adequate-precision arithmetic, adequate 
mesh densities (with FEM), and, where confidence or expe- 
rience is lacking, a thorough testing of the robustness of the 
solution against changes in the geometry or mesh density. 

Now, the energy stored in the resonator's electromagnetic 
field is U = (1/2) /// Ai|H| 2 dV, where H is the infinite-Q 
solution generated by the PDE-solver, and fi is the common 
permeability for the resonator's interior. For resonators that 
are axisymmetric, the 3D volume integral J J J dV over the 
resonator's interior reduces to the 2D integral J J (2-Kx)dxdy 
over its medial cross-section. The surface current induced in 
the resonator's enclosing perfect-electric wall is given by (see 
ref. [18], page 205, for example) J s = H t EnxH, where H t 
is the tangential component of H with respect to the resonator's 
electric-wall boundary. 

One now exploits (first-order) perturbation theory, and 
equates the current immediately stated above with that which 
would be induced into the electric walls of a resonator, 
identical to its loss-less antecedent, but for it having electric 
walls of finite conductivity. The equating of the two currents 
assumes that the lossy walls are nevertheless made out of 
(or coated with) a sufficiently 'good' conductor, such that the 
change from loss-less to lossy does not significantly affect the 
shapes of the resonator's modes. This will typically be the case 
for a cavity resonator exciting low-order modes at microwave 
frequencies, provided its walls are made from any standard 
(electrically good) metal, such as copper; again, see references 
[18] and/or [34] for further explanation/quantification. 

The time-averaged(-over-a-cycle) power lost by the res- 
onator through resistive heating in its imperfect electric walls 
is thereupon given by P loss = (1/2) J / i? s |n x H| 2 dS, where 
the 2D surface integral J J dS over the resonator's presumed 
axisymmetric enclosing boundary reduces to the ID integral 
J (27ra:)dl around the periphery of its medial (x-y) cross- 
section; R s = (nffi/a) 1 / 2 is the wall's surface resistivity, 
where a is the wall's (bulk) electrically conductivity, and / the 
mode's frequency. The quality factor, defined as 27r/ U/P\ oss , 
due to the wall's resistive losses can thus be expressed as: 

Qwaii = = (Anffxa^K, (26) 

where A, which has the dimensions of a length, is defined as 

a J7/|H| 2 dV 
J7|nxH| 2 dS 

JJ[(H*f + (H*r + (Hy) 2 ]xdxdy 
] x[\H<t>\ 2 + \Hvdn x - H x n y \ 2 ]&\ ' 1 ' 

Both integrals (numerator and denominator), hence Q wa ii 
itself, can be readily evaluated using the PDE-solver's post- 
processing features; explicit PDE-solver-ready forms of each 
integrand are provided in this paper's separate Appendix [37]. 
In the using of equation 26, it should be pointed out that, at 
liquid-helium temperatures, the bulk and surface resistances of 



EX-HOUSE 2D FINITE-ELEMENT SIMULATION OF WHISPERING-GALLERY MODES 



8 



metals can depend greatly on the levels of (magnetic) impu- 
rities within them [38], and the text-book /~ 1/2 dependence 
of surface resistance on frequency is often violated [39]. 

IV. Radiation Loss in Open Resonators 

A. Open resonators: preliminary remarks 

Many whispering-galley mode resonators (both microwave 
[40] and optical [12], [14]) of interest experimentally are 
not shielded by an enclosing metal wall: they are open. 
In consequence, the otherwise highly localized WG mode 
supported by the resonator spreads throughout free-space 7 , 
leading to the conveyance of energy away from the mode's 
bright spot (where the electric- and magnetic-field amplitudes 
are greatest) through radiation. Provided the equivalent closed 
resonator's enclosing electric (i.e. metal) wall is stationed 
sufficiently far out in the WG mode's evanescent zone, the WG 
mode's form in its near-field will be the same (to the degree 
of equivalence required here) in both the open-resonator and 
closed-resonator cases. One can thus calculate the mode's 
near-field form through the method developed in section II, 
or by some other method, as applied to the closed resonator; 
in particular, the electric and magnetic field strengths, E and 
H or, equivalently, the vector potential A, just outside of the 
surface(s) of the resonator's dielectric component(s) can be 
determined. Having done so, the WG mode's far-field form 
(i.e. its 'radiation pattern') in the case of the open resonator 
can be calculated by invoking the so-called Field Equivalence 
Principle [41], [42], where A or (the tangential components 
of) E and H over the said surface(s) are regarded, in Huygen's 
picture, as (secondary) sources radiating into free-space. The 
calculation can be implemented through a standard retarded- 
potential (Green function) approach [17], [42], incorporating 
(if necessary) a multipole expansion. The mode's radiative 
loss, hence Q, can be subsequently calculated from the radia- 
tion pattern determined by integrating Poynting's vector over 
all solid angles. With due care, the resulting estimate of the 
mode Q will be highly accurate. But such a program of work 
-for lack of novelty rather than utility- shall not be undertaken 
here. 

B. Estimators of radiation loss 

Instead, two different (but related) 'trick' methods for 
estimating the radiative Q of an open (dielectric) resonator 
are described here. As the first method underestimates the Q, 
while the second overestimates it, the two in conjunction can 
be used to bound the Q from below and above. Moreover, 
both can be implemented as straightforward 'add-ons' to the 
2D PDE-solver's computational environment, as already con- 
figured for solving closed loss-less resonators (as per section 
II). It should be added that these two methods are not restricted 
to axisymmetric resonators per se. 

7 As understood by Kippenberg [28], this observation implies that the 
support of equation 24's J J J h _ B ...dV integral, as it covers the WG mode's 
bright spot, must be somehow limited, spatially, or otherwise (asymptotically) 
rolled off, lest the integral diverge as the so-called 'quantization volume' 
associated with the radiation extends to infinity. 



1) Underestimator via (imperfect) retro-reflection: Con- 
sider an otherwise loss-less open resonator, supporting a 
spatially concentrated mode, i.e. one with a bright spot, 
that radiates into free-space. As stated above, the tangential 
electric and magnetic fields on any closed surface in the near- 
field surrounding this mode's bright spot can, by the Field 
Equivalence Principle, be regarded as a (secondary) source 
of this radiation. Now consider a closed, completely loss-less 
equivalent of the open resonator, formed by placing a cavity 
around it, whose enclosing perfect-electric wall lies in the said 
localized mode's radiation zone. It is noted here, for future 
reference, that this perfect electric wall will force the tangential 
component of the electric field strength to vanish everywhere 
on it, i.e. E — n(E • n) = 0, where n is the wall's unit normal 
vector. The above-mentioned secondary source generates an 
outward-going traveling wave wave which, but for the cavity, 
would lead to radiation. Instead, with the cavity in place, a 
standing (as opposed to traveling) wave arises. Now, suppose 
that the shape of the cavity's electric wall, and its location 
with respect to the source, is chosen to predominantly reflect 
the source's outward-going wave back to the source such that 
the resultant inward-going wave interferes constructively with 
the outgoing wave over the whole of the source's surface. In 
other words (ID analogy), on regarding the cavity as a short- 
circuit-terminated transmission line, one attempts to adjust the 
length of the line such that its input (analogous to the source's 
surface) is located at an antinode of the line's standing wave. If 
such a retro-reflecting (+ phase-length adjusted) cavity can be 
devised then, in particular, the measured/simulated tangential 
magnetic field, H t , just inside of the cavity's electric wall will 
be exactly twice that of the outward-going wave for the open 
resonator H[ ad ' at the same location -but without the cavity's 
electric wall in place. In practice, the source will not be located 
exactly at an antinode (over the whole of its surface) and 
thus H t > 2H£ . The radiative loss for the open resonator 
can be evaluated by integrating the cycle-averaged Poynting's 
vector corresponding to the outward-going wave's inferred 
tangential magnetic field over the electric wall's surface; i.e. 
-Prad. = (1/2) / / z |H[ ad | 2 dS, where z is the impedance of 
free space. A bound on the open resonator's radiative Q-factor 
can thus be expressed as 

Qrad. > (8tt// C )A, (28) 

where A is exactly that given by equation 27 but with the (loss- 
less) electric wall now in the radiation zone. Provided the PDE 
solver is able to accurately calculate the (faint) electromagnetic 
field on the rad.-zone cavity's enclosing boundary, it can 
again be readily determined. It is further remarked here that 
the above -admittedly rather heuristic and one-dimensional- 
argument, is strongly reminiscent of Schelkunoff 's induction 
theorem [41], [43], which is itself a corollary of the (already- 
mentioned) Field Equivalence Principle. Through analogy 
to this theorem, the author conjectures that equation 28 
holds equally well for fully vectorial waves (as governed by 
Maxwell's equations) in 3D space as for scalar waves along 
ID transmission lines -as argued above. 
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2) Overestimator via (imperfect) outward-going free-space 
impedance match: A complementary estimator to the one 
above can be constructed by replacing the above cavity's 
electric wall with an 'impedance-matched' one, where the 
tangential magnetic and tangential electric fields at every point 
on the wall are constrained so as to correspond to those of 
an outward-going plane traveling wave, propagating in the 
direction of the wall's local normal and in an outward-going 
sense. In other words (ID analogy), one attempts to confront 
the secondary source's outward going wave with a matched 
surface that reflects nothing back. For plane-wave radiation, 
this constraint can be expressed as z n X H = E n(E • n), 
where n is the surface's inward-pointing normal. Note that 
one does not constrain the direction (polarization) of E or 
H in the wall's local plane; one only demands that the two 
fields be orthogonal and that their relative amplitudes be in the 
ratio of the impedance of free space z . Upon differentiation 
w.r.t. time and using Maxwell's displacement-current equation, 
this relation can be re-expressed as 

V X H — n[(V X H) • n] — (l/c)n X — = 0. (29) 

For a given eigenmode, and generalizing somewhat, the con- 
straint can be recast as 

cos((9 mix ){V x H n[(V x H) • n]} 

+ sin(0 mix ) ic/nxH = 0, (30) 

where / is the mode's frequency (in Hz), c = 2tt/c as before, 
and # mix is a 'mixing angle'. In the impedance-matched case 
(outward going plane wave in free space), 6* m j x = 7r/4, and 
the above equation reduces to 

V x H - n[(V x H) • n] + i c/n x H = 0. (31) 

More generally, the first and second terms on the left-hand 
side of equation 30 can be viewed as representing electric- 
(cf. equation 3) and magnetic-wall (cf. equation 5) boundary 
conditions, respectively, where the latter corresponds to that 
used in subsubsection IV-B.l above. The boundary condition 
can be continuous adjusted between these two extrema by 
varying the mixing angle 6* m ; x ; for the sake of completeness, 
#mix = — 7r/4 corresponds to an inward-going, as opposed 
to an outward-going, impedance match. Note that, unless 
$mix = Ntt/2 for integer TV, the square root of minus one 
in equation 30 breaks the hermitian-ness of the matrix that 
the PDE-solver is required to eigensolve, leading to solved 
modes with complex eigenfrequencies ,/ mo dc- As exploited 
by Srinivasan et al [14], the inferred quality factor for such 
a mode due to radiation through/into its bounding wall is 
given by Q inf . = 3?[/ m odo]/23?[/ mo do], where &[...] and 9f[...] 
represent taking real and imaginary parts (of the complex 
eigenfrequency), respectively. 

Note that the accuracy of the method will depend on the 
degree to which the imposed surface impedance agrees with 
that of the true outward-going traveling wave, as generated 
by the open resonator's (secondary) source, over the chosen 
bounding surface. If the source were an infinitessimal(finite) 
multi-pole, then a surface in the form of a finite(infinite) sphere 
centered on the source, with the constraint 31 imposed on its 



surface, would perfectly match to the source's radiation (i.e. no 
traveling wave would get reflected back from it). In general, 
however, for a finite radiator, the chosen surface (necessarily 
of finite extent) will not lie everywhere normal to the outward- 
going wave's Poynting's vector and back reflections will 
result, leading to a smaller imaginary part in the simulated 
eigenmode's frequency, thus causing Qi„f. to overestimate the 
true radiative Q. Thus, one may state 

Qrad. < »[/mode]/29[/mode], (32) 

approaching equality on perfect matching. Again, the au- 
thor conjectures that, despite the rather heuristic and one- 
dimensional argument stated above, inequality 32 holds in 
general. Used together, equations 28 and 32 provide a bounded 
estimator on the true radiative Q. 

Comment: As alluded to at the beginning of this section, 
the author recognizes that alternative (one might argue rather 
more 'empirical') approaches, based on surrounding (cladding) 
the otherwise open resonator with sufficiently thick layers 
of impedance-matched absorber [i.e., with the absorber's di- 
electric constant set equal to that of free space except for 
a small imaginary part (loss tangent) causing the outward- 
going wave to be gently attenuated with little back reflection]. 
The 'boundary-alteration' method described above has the 
advantage of not extending the footprint of the PDE-solver's 
modeled region (thus not requiring the mesh of finite elements 
to be extended). 

V. Example applications 

The author has deployed the methodologies expounded in 
sections II through IV above to model several different sorts 
of resonator. Where possible, he chose resonators with shapes 
and properties that had already been published -so as to afford 
comparisons. Each of the characteristics considered in sections 
III and IV was evaluated for at least one such model resonator. 
The COMSOL applications (as '.MPH' files) that the author 
constructed to simulate these resonators can be obtained from 
him upon request. 

A. UWA 'sloping-shoulders' cryogenic sapphire microwave 
resonator 

This resonator, as designed and assembled by workers at the 
University of Western Australia (UWA) [44], [15], comprises a 
piece of monocrystalline sapphire mounted within a cylindrical 
metal can. The can's internal wall and the sapphire's outer 
surface exhibit rotational symmetry about a common axis. 
Furthermore, the optical (or 'c') axis of the sapphire crystal 
is, to good approximation, oriented parallel to this geometric 
axis. The resonator can thus be taken (and modeled) as being 
electromagnetically axisymmetric. The sapphire piece's medial 
cross-section (one half thereof) is shown in Fig. 2(a). What 
makes the resonator awkward to simulate accurately via the 
semi-analytic MM-SV method [15] is its sloping shoulders (Si 
and S2 ibid.), whose surface normals are neither purely axial 
nor purely radial. 
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a: 



c: 




Fig. 2. UWA 'sloping-shoulders' cryogenic sapphire resonator: (a) me- 
dial cross-section through it; the grey(white) shading corresponds to sap- 
phire(vacuum); Si and S2 indicate the sapphire piece's two (upper and lower) 
'shoulders', (b) the FEM-based PDE-solver's meshing of the resonator's 
model structure; for clarity only every other meshing line is drawn [i.e. (b) 
displays the 'half-mesh']; (c) magnetic field intensity of the resonator's 
WG£i4 I o,0 mode; (d) electric field intensity of the same; in both (c) and (d) 
faint white arrows indicate the direction of the magnetic and electric field, 
respectively, in the medial plane; 'intensity' here means the absolute value 
of the vectorial H (or, equivalently, B) field, displayed on a logarithmic grey 
scale -darker being more intense. 



The resonator's form, as the author encoded it into the PDE- 
solver, is based on figure 3 of ref. [15] 8 . For the simulation 
presented here, the author took the piece's outer diameter, the 
length of its outer axial sidewall, the axial extent of each 
sloping shoulder, and the radius of each of its two spindles 
to be, at liquid-helium temperature {i.e. the dimensions here 
stated include cryogenic shrinkages -see section VI) 49.97, 
19.986, 4.996, and 19.988 mm, respectively. The sapphire 
crystal's cryogenic permittivities were taken to be {ej_, £||} = 
{9.2725,11.3486}, as stated in ref. [22]. Note that, though 
coaxial, the sapphire piece and the can do not exactly share a 
common transverse ('equatorial') mirror plane, thus precluding 
any speeding up of the simulation through the placement (in 
the model) of a magnetic or an electric wall on the equatorial 
plane, thereupon halving the 2D region to be analyzed 9 . 

Fig. 2(b) displays the FEM-based PDE-solver's meshing of 
the resonator structure. Here, the resonator's interior dielec- 
tric domains were regularly decomposed into quadrilaterals 
(as opposed to triangles), with no quadrilaterals straddling 
interfaces between different materials. The mesh comprised, 

8 It is commented parenthetically here that the shape of the sapphire piece 
in figure 3 of ref. [15] is not consistent with the dimensions stated in the 
same: its outer axial sidewall is too long and the slope of its shoulders too 
small with respect to the stated axial dimensions. 

'These two boundary conditions would lead to symmetric (N) and anti- 
symmetric (S) solutions, respectively -see ref. [45]. 



TABLE I 

ELECTRIC FILLING FACTORS FOR THE WGEl4,o,0 MODE OF THE UWA 







RESONATOR 




pdir. 


radial 


azimuthal 


axial 


sapphire 


0.80922 


0.16494 


7.016 X 10~ a 


vacuum 


0.01061 


8.0533 x 10~ 3 


1.6543 x 10~ 4 



in COMSOL's vernacular 10 , 7296 base-mesh elements and 
88587 degrees of freedom ('DOF'). It typically took around 
75 seconds, to obtain the resonator's lowest (in frequency) 16 
modal solutions, for a single, given azimuthal mode order M, 
at [with respect to Fig. 2(b)] full mesh density, on a standard, 
2004-vintage personal computer (2.4 GHz, Intel Xeo CPU), 
without altering the PDE-solver's default eigensolver settings. 
With the azimuthal mode order set at M = 14, the model 
resonator's WGEi4, ,o mode was found to lie at 11.925 GHz, 
to be compared with 11.932 GHz found experimentally [15]. 
Wall loss: This mode's characteristic length A was determined 
to be 2.6 x 10 4 . Based on ref. [39], one estimates the surface 
resistance of copper at liquid-helium temperature to be 7 x 
10 -3 per square at 11.9 GHz, leading to a wall-loss Q of 
3.5 x 10 11 for the WGEi 4 , ,o mode. As this is at least an order 
of magnitude greater than what is observed experimentally, 
one concludes that wall losses do not substantially limit the 
UWA resonator's experimental Q. 

Filling factor: Using equation 25, the electric filling factors 
for the WGE14 0.0 mode have been evaluated. The results, 
presented in TABLE I, are to be compared with those stated 
in Table IV of ref. [15]: they are in good agreement with 
those loc. cit. (labeled 'FE'), which were obtained via a wholly 
independent finite-element analysis. 

B. Toroidal silica optical resonator [Caltech] 

The resonator modeled here, based on ref. [13], comprises a 
silica toroid, where this toroid is supported through an integral 
interior 'web', such that the toroid is otherwise suspended in 
free space above the resonator's substrate. This arrangement is 
shown in Fig. 3(a). The toroid's principal and minor diameters 
are set at {D,d} = {16,3} pm, respectively. The silica di- 
electric is presumed to be wholly isotropic (i.e., no significant 
residual stress) with a relative permittivity of e s n. = 2.090, 
corresponding to a refractive index of n s i\. = ^/e s n. = 1.4457 
at the resonator's operating wavelength (around 852 nm) and 
temperature. The FEM model's mesh covered an 8-by-8 /im 
square [shown in dashed outlined on the right of Fig. 3(a)] 
in the medial half-plane containing the silica toroid's circular 
cross-section. A pseudo-random triangular mesh was gener- 
ated (automatically) with an enhanced meshing density over 
the silica circle and its immediately surrounding free-space; 
in total, the mesh comprised 5990 (base-mesh) elements, 
with DOF = 36279. Temporarily adopting Spillane et aVs 

10 The size/complexity of a finite-element mesh is quantified, within COM- 
SOL Multiphysics, by (i) the number of elements that go to compose its 
so-called 'base mesh' and (ii) its total number of degrees of freedom ('DOF') 
-as associated with its so-called 'extended mesh'. Consult the package's 
documentation for further clarifications. 
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terminology, the resonator's fundamental TE-polarized 93rd- 
azimuthal-mode-order mode (where by 'TE' it is here meant 
that the polarization of the mode's electric field is predomi- 
nantly aligned with the toroid's principal axis -not transverse 
to it) was found to have a frequency of 3.532667 x 10 14 Hz, 
corresponding to a free-space wavelength of A = 848.629 nm 
(thus close to 852 nm). Using this paper's equation 24, this 
mode's volume was evaluated to be 34.587 /jm 3 ; if, instead, 
the definition stated in equation 5 of ref. [13] is used, the 
volume becomes 72.288 /im 3 —i.e. a factor of greater. 
These two values straddle (neither agree with) the volume 
of approx. 55 /im 3 , for the same dimensions of silica toroid 
and the same (TE) mode-polarization, as inferred by eye-and- 
ruler from figure 4 of ref. [13]. The author cannot explain the 
discrepancy. 

It is pointed out here that the white arrows in Fig. 3 (at 
least those not anchored on the equatorial plane) are slightly 
but noticeably oriented away from vertical, indicating that 
the orientation of the mode's (vectorial) electric field is not 
perfectly axial -as per the transverse approximation taken in 
references [15], [28]. In other words, the arrows' lack of 
verticality reveals the inexact- or quasi-ness of the mode's 
transverse-electric nature, despite the mode's relatively high 
azimuthal mode order (I = M = 93). 



The model 'microdisk' resonator analyzed here, as depicted 



a: 




D 








Fig. 3. (a) Geometry (medial cross-section) and dimensions of a model 
toroidal silica microcavity resonator -after ref. [13]; the torus' principal 
diameter D = 16 fim and its minor diameter d = 3 /im; the central 
vertical dashed line indicates the resonator's axis of continuous rotational 
symmetry, (b) False-color surface plot of the (logarithmic) electric-field 
intensity |E| 2 within the dashed box appearing in (a) for this resonator's 



TE 



p— 1 ,m 



=93 whispering-gallery mode. White arrows indicate the electric 



field E's magnitude and direction in the medial plane. 



C. Conical microdisk optical resonator [Caltech ] 

The mode volume can be reduced by going to smaller 
resonators, which, unless the optical wavelength can be com- 
mensurately reduced, implies lower azimuthal mode order. 



a: 





Fig. 4. (a) Geometry (medial cross-section) and (half-)meshing of model 
GaAlAs microdisk resonator -after ref. [14]; the disk's median diameter is D 
= 2.12 fim and its thickness (axial height) t = 255 nm; its conical external 
sidewall subtends an angle 8 = 26° to the disk's (vertical) axis; for clarity, 
only every other line of the true (full) mesh is drawn. The modeled domain 
in the medial half-plane is a rectangular stretching from 0.02 to 1.5 (im in the 
radial direction and from -5 to +5 (im in the axial direction, (b) False-color 
surface plot of the (logarithmic) electric-field intensity |E| 2 for the resonator's 
TE p= i. m= ii mode at A = 1263.6 nm. Again, white arrows indicate the 
electric field's magnitude and direction in the medial plane. 

in Fig. 4(a) is the author's attempt at duplicating the structure 
defined in figure 1(a) of Srinivasan et al [14]; as in their 
model, the disk's constituent dielectric (in reality, layers of 
GaAs and GaAlAs) is approximated as a spatially uniform, 
isotropic dielectric, with a refractive index equal to n = 
3.36. The FEM-modeled domain in the medial half-plane 
was divided into 4928 quadrilateral base-mesh elements, with 
DOF = 60003. Adopting the same authors' notation, the 
resonator's TE p= i. m= n whispering-gallery mode, as shown 
in Fig. 4(b), was found at 2.372517 x 10 14 Hz, equating to 
A = 1263.6 nm; for comparison, Srinivasan et al found an 
equivalent mode at 1265.41 nm [as depicted in their figure 
1(b)]. The white arrows' lack of verticality in this article's 
Fig. 4(b) implies that the orientation (i.e. polarization) of 
the magnetic field associated with the true, quasi-TE p= i m= n 
mode deviates significantly from axial (as would be the case 
within a transverse approximation). 

Mode volume: Using this paper's equation 24, but with the 
mode excited as a standing-wave (doubling the numerator 
while quadrupling the denominator), the mode volume is 
determined to be 0.1484 x ^m 3 ~ 2.79(A/n) 3 , still in good 
agreement with Srinivasan et al's ^ 2.8(A/ra) 3 . 
Radiation loss: The TE p= i )m= n mode's radiation loss was es- 
timated by implementing both the upper- and lower-bounding 
estimators described in subsection IV-B. Here, the microdisk 
and its mode were modeled within a near- spherical volume 
(equating to a half-disk in the medial half-plane, with a 
semicircle for its outer perimeter), on whose outer boundary 
different electromagnetic conditions were imposed -see Fig. 5. 
With an electric-wall condition (i.e. equations 2 and 3 or, 
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Fig. 5. Radiation associated with the same [TE p= i im= n, A = 1263.6 nm] 
whispering-gallery mode as presented in Fig. 4; here, false-color maps of the 
squared magnitude of the mode's magnetic field strength are plotted out to 
the modeled domain's near-spherical outer boundary, corresponding in the 
medial half-plane to a semi-circle 12 fim in diameter, whose center lies at 
a radial coordinate of 0.01 /im on the microdisk's mid-plane. [In reality, the 
microdisk's substrate would occupy a considerable part of the meshed half- 
disk's lower quadrant, but the model here assumes that, with the exception 
of the microdisk itself (a dielectric), both quadrants are filled with free space 
—into which the the whispering gallery mode radiates.] All three maps use 
the same absolute false-color scale, (a) standing-wave (equal outward- and 
inward-going) radiation with the outer semicircular boundary set as a magnetic 
wall; (b) the same but now with the boundary set as an electric-wall; (c) 
somewhat traveling (more outward- than inward-going) radiation with the 
boundary's impedance set to that of an outward-going plane-wave in free 
space (and with the normal magnetic field constrained to vanish). That (c)'s 
radiation field is somewhat dimmer than (b)'s is consistent with the different 
estimates of the resonator's radiative Q corresponding to (a)-(c) [see text]. 



equivalently, 17 and {18, 20}) imposed on the volume's whole 
boundary [as per Fig. 5(b)], the right-hand of equation 28 
was evaluated. And, with the E x n = condition (viz. 
equation 3) on its outer semi-circle replaced by the outward- 
going-plane-wave(-in-free-space) impedance-matching condi- 
tion (viz. equation 30), while the H • n = condition (equation 
2) is maintained, the right-hand side of equation 32 was 
evaluated for the radiation pattern displayed in Fig. 5(c). For a 
pseudo-random triangulation mesh comprising 4104 elements, 
with a DOF of 24927, the PDE solver took, on the author's 
office computer, 6.55 and 13.05 seconds, corresponding to 
Figs. 5(b) and (c), respectively 11 , to calculate 10 eigenmodes 
around 2.373 x 10 14 Hz, of which the TE p=1 m=11 mode 
was one. Together, the resultant estimate on the TE p=1 m=11 
mode's radiative-loss quality factor is (1.31 < Q ra d. < 
3.82) x 10 7 , to be compared with the estimate of 9.8 x 10 6 (at 
1265 nm) reported in table 1 of ref. [14] 12 . The standing-wave 
radiation field in Fig. 5(b) could have been made dimmer (thus 
increasing A, hence the inferred Q) by adjusting ('tuning') the 
meshed half-disk's diameter -so as to put the microdisk/near- 
field TE p= i im= n mode, viewed as a secondary source of 
radiation, closer to an antinode of the cavity's standing- 
wave field. Also, the simulations associated with Fig. 5 could 

1 1 The complex arithmetic associated with the impedance-matching bound- 
ary condition meant that the PDE solver's eigen-solution took approximately 
twice as long to run with this condition imposed -as compared to the 
electric- (or magnetic-) wall boundary conditions that do not involve complex 
arithmetic. 

12 The author chose the diameter of the outer semicircular boundary in 
Figs. 5(a)-(c) arbitrarily to be 12.0 fim in advance of knowing what upper 
and lower bounds on Q ra( j. such a choice would give; he did not subsequently 
adjust the diameter and/or shape of this boundary to bring the bounds any 
closer together. 



certainly have run (with tolerable execution times) on a denser 
finite-element mesh. 

D. 3rd-order Bragg-cavity alumina:air microwave resonator 

Commercial FEM-based PDE-solvers (viz. the COM- 
SOL/FEMLAB package used by the author for this article) 
permit the simulation of arbitrarily complex structures and, 
moreover, provide efficient languages and tools for represent- 
ing and constructing (and modifying) them. Through such 
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Fig. 6. (a) Geometry (medial cross-section) of a alumina:air 3rd-order Bragg- 
cavity resonator within a cylindrical metallic can (electric walls); the can's 
interior surfaces are represented by a solid black line; its interior diameter 
equals its interior height (and thus this black line takes the form of a square); 
the horizontal and vertical grey (or pink -in color reproduction) stripes 
denote cylindrical plates and barrels, respectively, of alumina; white squares 
correspond to regions of free-space (either air or vacuum); the vertical arrow 
indicates the resonator's axis of rotational symmetry; the dashed horizontal 
line (cf. Ml in Fig. 1) denotes a plane of mirror symmetry, on which an 
electric or magnetic wall is imposed, (b) False-color plot of the (logarithmic) 
electric-field intensity |E| 2 for a zeroth-azimuthal-mode-order (M = 0) mode 
at 8.0873 GHz, localized towards the resonator's center (bottom left in figure); 
(b) the same but for a sixth-azimuthal-mode-order (M = 6) mode at 20.0267 
GHz, strongly localized in the radial directions but less so in the vertical 
direction. 

a PDE-solver, the method described in sections II through 
IV can be applied to axisymmetric dielectric resonators of 



BX-HOUSE 2D FINITE-ELEMENT SIMULATION OF WHISPERING-GALLERY MODES 



13 



arbitrarily complex medial cross-section, the only require- 
ments being that each such cross-section is (i) bounded (either 
externally by an enclosure or internally as an excluded region, 
or both) by metallic walls and (ii) decomposable into definable 
regions of uniform dielectric. This ability to cope with struc- 
tural complexity is exemplified here in a modest way through 
the simulation of a 3rd-order Bragg-cavity alumina: air mi- 
crowave resonator whose geometry is shown in Fig. 6(a). This 
resonator's model geometry was generated straightforwardly 
through a script written in MATLAB. The resultant FEM 
mesh in COMSOL comprised 4356 base-mesh elements, with 
53067 degrees of freedom (DOF), corresponding to 12 edge 
vertices per A/4 interval of air [i.e., across each white square 
in Fig. 6(a)], and 6 vertices per A/4 interval of alumina [i.e., 
across each grey/pink 'strip', ibid.]. Figs. 6(b) and (c) display 
two different calculated modes that this resonator supports. 

VI. Determination of the permittivities of 

CRYOGENIC SAPPHIRE 

The author harnessed the method of simulation constructed 
in sections II and III to extract an independent determina- 
tion of the two dielectric constants of pure (HEMEX [46]) 
monocrystalline sapphire at liquid-helium temperature, based 
on some existing experimental data [47]. This data, as is listed 
in the four right-most columns of TABLE II, comprised 13 : the 
centre frequencies, FWHM widths, turnover temperatures, and 
'Kramers' splittings for a set of 16 resonances, as measured 
on a (one of a pair of) cryogenic sapphire resonator(s), as 
shown, without its enclosing can, in Fig. 7(a). Only two 
resonances out of this set (viz. Nln and S2g) had hitherto 
been identified -via MAFIA [7], [9] simulations 14 . This cryo- 
sapphire resonator's complete, detailed model geometry, as 
shown in Fig. 7(b), was coded into a MATLAB script. This 
script contained, for example, the dimensional parameters 
specifying the form of the sapphire ring's (large) external 
and (smaller) internal chamfers. The model geometry took 
into account the shrinkages of the resonator's constituent 
materials from room temperature (293 K) down to liquid- 
helium temperature (4.2 K). The two cryo-shrinkages of 
sapphire (a uniaxial crystal) were calculated by integrating 
up [48] the linear-thermal-expansion data stated in Table 4 of 
ref. [49] (identical to that stated in TABLE 1 of ref. [50]): 
(1.0 - 7.21 x 10~ 4 ) and (1.0 - 5.99 x 10~ 4 ) for directions 
parallel and perpendicular to sapphire's c-axis, respectively. 
The cryo-shrinkage of (isotropic) copper was taken directly 
from Table F at the back of ref. [51] 15 : (1.0 - 3.26 x 10~ 3 ). 
The values of sapphire's two dielectric constants were initially 
set equal to those specified in ref. [22]: e± = 9.2725 and 
en = 11.3486 16 . Fig. 7(b)'s geometry was meshed with 

13 Though not listed in TABLE II, the measured insertion loss (i.e. S21 at 
line center) for each resonance was also available. 

14 The as-measured resonator was developed as part of a local 'flywheel' 
oscillator for supplying NPL's Cs-fountain(s) with an ultra-frequency-stable 
9.1926 ... GHz reference, with the resonator operating on (as it turned out) 
the S2g WG mode. 

15 Ref. [52] provides linear-thermal-expansion data for copper as a function 
of temperature -useful for design purposes. 

16 These values are consistent with e±_ = 9.272 and ey = 11.349, as stated 
in ref. [15]. 




Fig. 7. (a) Close-up of one of NPL's two (nominally identical) Cs-fountain 
cryo-sapphire resonators, with its outer copper can removed. The resonator's 
chamfered HEMEX sapphire ring has an outer diameter of ~46.0 mm and a 
axial height of ~25. 1 mm. This ring's integral interior 'web', 3mm thick, is 
oriented parallel to and centered (axially) on the ring's equatorial plane; the 
web is supported through a central copper post, which is in turn connected 
[indirectly -through a thin, annular stainless steel 'shim' (not visible)] to 
the resonator's copper lid (onto which the removed can is secured). Note 
that the sapphire's high refractive index falsely exaggerates [cf. the true 
relative dimensions shown in (b)] the ring's internal diameter of ~20.0 mm. 
Above the ring lie two loop probes for coupling, electromagnetically, to 
the resonator's operational whispering-gallery mode. [As finally configured, 
these probes were withdrawn upwards several mm's closer to the lid and 
thus further (axially) from the ring], (b) geometry of the resonator in medial 
cross-section; pink/grey indicates sapphire, white free space; bounding these 
dielectric domains, and shown as thick solids lines, are copper surfaces 
belonging to the resonator's can, lid and ring-supporting post; the resonator's 
cylindrical axis (r or x = 0) is shown as a dashed vertical line, (c) false- 
color map (logarithmic scale) of the magnetic (H) field's squared magnitude 
for the resonator's 1 lth-azimuthal-mode-order fundamental quasi-transverse- 
magnetic (Nln in ref. [45]'s notation) whispering-gallery mode at 9.146177 
GHz (simulated), as detailed on the 6th row of TABLE II. The white arrows 
indicate the magnitude and direction of this mode's electric (E) field in the 
medial plane. 



quadrilaterals over the medial half-plane, with 8944 elements 
in its base mesh, and with DOF = 108555. [These quadri- 
laterals followed the sloping chamfers of the resonator by 
taking the shape of trapezoids.] For a given azimuthal mode 
order M, the calculation of the lowest 16 eigenmodes took 
around 3 minutes on the author's office PC (as previously 
specified). Fig. 7(c) shows the form of the resonator's Nln 
whispering-gallery mode, corresponding to the 6th row of 
TABLE II. Filling factors were then calculated to quantify 
each frequency's sensitivity to changes in the sapphire's two 
dielectric constants (en and e±). The author identified each 
of the 16 experimental resonances with a particular simulated 
WG mode, aiming to minimize the residual (simulated-minus- 
measured) frequency difference (i.e. the sum \ 2 variance over 
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TABLE II 

NPL'S CRYOGENIC SAPPHIRE RESONATOR: SIMULATED AND 
EXPERIMENTAL WG MODES COMPARED 
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"the nomenclature of ref. [45] is used for this column. 
' J full width half maximum (-3 dB) 

c the difference in frequency between the orthogonal pair of standing- 
wave resonances (somewhat akin to a 'Kramers doublet' in atomic physics) 
associated with each WG model; the experimental parameters stated in other 
columns correspond to the strongest resonance (greatest S21 at line center) 
of the pair. 



the left-most column in TABLE II), whilst requiring that the 
other measured attributes (e.g. insertion loss, linewidth) of 
the resonances identified to the same 'family' of WG modes 
(e.g. SI or N2) varied smoothly with the azimuthal mode 
order M . With the identifications of the experimental modes 
'locked' as per the 4th column of TABLE II, the model's two 
sapphire dielectric constants were adjusted from their initial 
values to minimize x 2 ('least squares'). The resultant best-fit 
values were: 



e ± = 9.285 (±0.010); 
6,1 = 11.366 (±0.010). 



(33) 
(34) 



With the two dielectric constants set to these values, the WG 
modes were recalculated (an in-principle superfluous check); 
the first three columns in TABLE II result from this recalcu- 
lation (the filling factors hardly changed from their original 
fitted values). Pending the construction of a more detailed 
and precise error budget, the provisional ±0.010 uncertainty 
assigned to the values of both dielectric constants reflects 
their observed shifts upon refitting with a few 'problematic' 
experimental modes identified with different simulated ones 11 . 

I7 A few of the experimentally measured modes had a number simulated 
modes in close proximity to their center frequencies; note that the centers 
of several circles lie close to certain horizontal lines in see Fig. 8. After 
considering all available pieces of experimental information (viz. linewidth, 
insertion loss, and turnover temperature), doubt still remained as to the correct 
identifications for some of them. The 4th column of TABLE II represents the 
most likely, but not the only conceivable, set. Though a 'prettier' (and perhaps, 
even, more accurate) determination could have been presented by dropping 
these problematic modes/identifications from the least-squares fit, the author - 
given the purposes of this paper- elected to fit all 16 experimentally measured 
resonance, keeping the generic problem of mode identification to the fore.] 



Compared to these identification-related shifts, the systematic 
errors associated with a finite meshing density -as analyzed 
quantitatively in ref [33]), or the experimental uncertainties 
associated with the resonator's geometric shape (particularly 
the diameter and height of the sapphire ring) were quite 
negligible. The specified (i.e. contracted) tolerance on the 
alignment of the sapphire crystal's c-axis with respect to the 
geometric axis of the 'cored' cylinder from which the two 
rings were cut (through orientation-preserving methods) was 
only < 0.5 degrees. Though the effects of crystal misalignment 
cannot be modeled quantitatively with the method presented in 
this paper, for which axial symmetry is a requirement, it can be 
estimated that, given the ^two-parts-in-ten contrast between 
the sapphire's parallel and transverse permittivities, such a 
misalignment should make a significant/dominant contribution 
to the 1-part-in-a-thousand residuals between the simulated 
and measured center frequencies (and thus the determination 
of and ey.) Even with azimuthal mode orders of M ~ 10, 
as is the case here, the narrowness of the measured WG 
resonances' Kramers splittings (listed in the right-most column 
of TABLE II, and generally less than 1 part per million 
relative to the absolute frequency) would indicate a much 
higher degree of rotational invariance, however. Though noting 
that center-frequency residuals of a few parts per thousand are 
not untypical for FEM-based simulations of WG modes [33], 
the author has yet to reconcile, convincingly, the residuals 
with their cause(s) -as would be required to construct a more 
detailed error budget. 

VII. Conclusion 

This paper demonstrates, through the explicit statement 
of weak-form expressions and boundary constraints, how a 
commercial (FEM-based) PDE-solver can be configured to 
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Fig. 8. Plot used to identify experimentally measured with simulated WG 
modes. Solid horizontal lines (16 in total) indicate the center frequencies of 
the former. Solid circles indicate the identification of a simulated mode with an 
experimental one (the difference in their frequencies corresponds to much less 
than a circle's radius in all cases); hollow circles indicate simulated modes that 
were not identified with any experimentally measured one. Quasi-transverse- 
magnetic (q-TM) and quasi-transverse-electric (q-TE) WG modes of the same 
family are joined by (blue-)dashed and (red-)dotted lines respectively; a few 
of the lowest-lying mode families are labeled using standard notation [45]. 



EX HOUSE 2D FINITE-ELEMENT SIMULATION OF WHISPERING-GALLERY MODES 



15 



simulate, quickly and to high accuracy, the whispering-gallery 
modes of axisymmetric dielectric resonators on standard com- 
puter hardware. The source codes/configuration scripts used 
to implement the simulations presented in section V of this 
paper are freely available from the author. 
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Appendix I 

Configuration of COMSOL Multiphysics for 

SIMULATING AXISYMMETRIC DIELECTRIC RESONATORS: 
EXPLICIT WEAK-FORM EXPRESSIONS 

It is here explained, in some detail, how to set up a 
dielectric-resonator simulation in COMSOL Multiphysics [10] 
-from scratch. These explanations should also be helpful 
to anyone wishing to modify one of the author's existing 
models -as incarnated in an .MPH file. At least in the first 
instance, it is recommended that the following instructions be 
meticulously adhered to -lest one stray from a tried-and-tested 
path. And it is suggested that the reader work through them 
with COMSOL Multiphysics open and running on his/her 
desktop. All menu items, expression names and variables 
associated with the program are displayed in typed text 
font. A good deal of supplementary information can be found 
in the documentation supplied with COMSOL Multiphysics 
itself; the author found the following chapters therein to be 
the most useful/relevant: 'PDE Modes for Equation-Based 
Modeling', 'The Weak Form', and 'COMSOL Multiphysics 
Scripting'. Upon reading these chapters, one might be left 
with the impression that COMSOL is simply not sufficiently 
flexible to embrace the task in hand (i.e. to implement sec- 
tions II through IV of this article explicitly); the following 
instructions demonstrate how COMSOL Multiphysics can, 
despite these first impressions, and most straight-forwardly, 
be so configured to implement the 2D simulation of isotropic 
dielectric resonators. From the beginning then: 

A. Setting up -fundamentals 

Get COMSOL Multiphysics up and running. Access the 
Model Navigator panel via File => New . . . and 
select the New tab if not already selected. 

(a) Select '2D' from the Space dimension: drop-down 
menu [note: do not choose 'Axial symmetric (2D) ']. 

(b) Browse to and select 'COMSOL Multiphysics 
=> PDE Modes => Weak Form, Subdomain' from the 
Application Mode navigator. 

(c) Type (verbatim) 'Hrad Hazi Haxi' into the 
Dependent variables : text field. These three variables 
are the radial, azimuthal and axial components of the magnetic 
field strength, respectively; all three are dependent on (i.e. are 
functions of) the Cartesian coordinates for the COMSOL 
simulation's 2D space, namely x (horizontal on the screen) 
and y (vertical) -both in units of metres [m]. The coordinate 
names 'x' and 'y' are already fixed by COMSOL (i.e. they 
are reserved symbols) and need not be explicitly entered (in 
COMSOL terminology, x and y are 'geometric variables'). 

(d) For the Application mode name: (default u) one 
can type in anything one likes. 

(e) Select 'Lagrange - Quadratic' from the 
Element : drop-down menu. [This choice is proven to 
work.] 

B. Constants 

All of the various constants (i.e. independent of x or 
y) included within the weak-form expressions given below 



are defined and described in TABLE III. The equivalent of 
this table needs to be typed (or loaded) into COMSOL's 
Options => Constants .... Each Expression thus 
Value therein [except those for eO, eperpO, eparaO - 
which define the (unit) relative permittivity of free-space], can 
be user-varied. But every Name should be entered verbatim; 
i.e., each constant must be named exactly as it appears in the 
expressions that subsequently include it. 

C. Expressions (for Postprocessing) 

The post-processing of the calculated magnetic-field 
strength (as a function of position) for each solved eigen- 
function is facilitated through the various definitions presented 
here. 

1) Scalar expressions: The equivalent of TABLE IV 
(or some subset thereof) needs to be typed into 
COMSOL's Options => Expression ... =4> 
Scalar expressions .... 

2) Subdomain expressions: The functionality of 
Subdomain expressions is required for generating 
post-processed fields, like the electric field strength E -as per 
the 6th, 7th and 8th entries in TABLE IV. Those constants 
associated with each such field's definitions, like (in the case 
of E) the relative permittivities epara and eperp, vary 
from one subdomain within the medial half place to another. 
The variation of these subdomain-dependent 'constants' 
is represented through Options => Expressions => 
Subdomain expressions; therein, the Name of each 
such variable is the same in each and every Subdomain 
(as identified by an integer), but its Expression reflects 
the variable value in the selected Subdomain. Thus, the 
Expression for epara in a Subdomain corresponding 
to (cryogenic and axisymmetrically oriented) sapphire would 
be eparal, with eparal defined (globally) as 11.3486 (or 
whatever) through TABLE III, whereas in a Subdomain 
corresponding to free space, the Expression for epara 
should be set to 1. Similarly (and more simply), the single 
Subdomain-dependent variable erel can be used to 
represent the variation of relative permittivity within an 
axisymmetric resonator containing solely isotropic dielectrics 
(incl. free space). 

D. Weak-form expressions 

The simulation's defining weak-form expressions are set 
up through the Physics =^> Subdomain Settings . . . 
control panel. On the left of this panel, first select the Groups 
tab. A New Group must be named and defined for each dielec- 
tric within the resonator being simulated. The author chose to 
name these dielectric Groups 'dielectric_0 : vacuum', 
'dielectric_l', 'dielectric_n', ... . For each 
dielectric Group, (in general) dielectric_n say, cor- 
responding weak-form expressions need to be entered into 
the weak terms (i.e. three slots or text fields), for ex- 
pressions involving spatial derivatives, and also into the 
dweak terms, for expressions involving temporal deriva- 
tives; these slots are accessed through the weak and 
dweak tabs, respectively, located on the right of the 
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Subdomain Settings . . . control panel. These terms 
govern the electromagnetic field in regions filled with the 
n-th dielectric, as specified in the Constants ... ta- 
ble introduced above. No other fields on the right of 
Subdomain Settings . . . need(/should) be touched; in 
particularly, don't monkey with the contr-tabbed sub-panel. 
Note that it is imperative that the Name of each constant 
entered into Options => Constants . . . above match 
(verbatim) its appearances within the expressions (below) that 
are entered into the weak and dweak text fields here. For 
each dielectric Group, two weak and one dweak terms are 
required: (i) a 'Laplacian' term (corresponding to the left most 
term on the left-hand side of equation 1 and (ii) a 'penalty' 
term, included to suppress spurious modes, corresponding to 
the middle term of the same. 

1) Laplacian term [first weak-term slot]: : The form of the 
Laplacian weak term, (V X H )- ( V x H), here given for the 
1st axisymmetric dielectric, is 18 



( (eperpl* (test (Hazi) *Hazi 

-M* (test (Hazi) *Hrad 

+Hazi*test (Hrad) ) 

+M~2*test (Hrad) *Hrad) 

+eparal *M" 2 *test (Haxi) *Haxi) /x 

+eperpl* (test (Hazix) * (Hazi-M*Hrad) 

+Hazix* (test (Hazi) -M*test (Hrad) ) ) 

-eparal*M* (test (Haxi) *Haziy 

+Haxi*test (Haziy) ) 

+x* (eperpl*test (Hazix) *Hazix 

+eparal* ( (test (Haxix) 

-test (Hrady) ) * (Haxix-Hrady ) 

+Haziy*test (Haziy) ) ) 

) / (eparal *eperpl ) 



(35) 



where Hazix denotes the partial derivative of Hazix with 
respect to the coordinate x, Hrady the partial derivative 
of Hazi with respect to y, etc.; test (Hazi) denotes 
the 'test function' of Hazi, etc. Its equivalent for the 2nd 
axisymetric dielectric is obtained by replacing eperpl by 
eperp2 and eparal by epara2, and so forth for all other 
axisymetric dielectrics (should more be required). The above 
expression can be significantly simplified for the (subdomain) 
Groups corresponding to isotropic dielectrics or free space 
(viz. dielectric_0); for computational efficiency, it is rec- 
ommended that these simplifications be implemented wherever 
possible. The required form of the Laplacian weak term, 
[(V X H ) • (V X H)]/ei, for the 1st isotropic dielectric is 



given explicitly as 

( (test (Hazi) *Hazi 

-M* (test (Hazi) *Hrad 

+Hazi*test (Hrad) ) 

+M"2* (test (Hrad) *Hrad 

+test (Haxi) *Haxi) ) /x 

+ (test (Hazix) * (Hazi-M*Hrad) 

+Hazix* (test (Hazi) -M*test (Hrad) ) ) (36) 

-M* (test (Haxi) *Haziy 

+Haxi*test (Haziy) ) 

+x* (test (Hazix) *Hazix 

+ ( (test (Haxix) 

-test (Hrady) ) * (Haxix-Hrady) 

+Haziy*test (Haziy) ) ) ) /el , 

where el is the material's dielectric constant (as appearing in 
TABLE III). The Laplacian weak term for the vacuum is 
the same with ei — > 1, and those for other isotropic dielectrics 
are similarly obtained by swopping ei with €2, £3, and so forth. 

2) Penalty (divergence-suppressing) term [second weak- 
term slot]: : The form of the penalty weak term, a(V • 
H ) • (V • H), the same for each subdomain Group, is 

alpha* ( (test (Hrad) *Hrad 
-M* (test (Hazi) *Hrad 
+Hazi*test (Hrad) ) 
+M~2*test (Hazi) *Hazi) /x 
+ (test (Haxiy ) 

+test (Hradx) ) * (Hrad-M*Hazi) 
+ (test (Hrad) -M*test (Hazi) ) 
* (Hradx+Haxiy ) 
+x* (test (Hradx) 
+test (Haxiy) ) * (Hradx+Haxiy) ) 

here, the coefficient alpha (whose value is determined 
through COMSOL's equivalent of TABLE III) controls the 
aggressiveness of the divergence suppression induced by this 
term. The remaining, 3rd slot, should be zero-filled. [As a 
general rule, unused weak-form slots should always be filled 
with zeroes -this applies to the the dweak term slots below.] 

3) Frequency term [first dweak-term slot]: The form of the 
temporal-derivative/frequency (so-called 'dweak') term H • 
d 2 H/d 2 t, common to all subdomain Groups, is entered into 
the first slot within the dweak-tabbed panel of Physics => 
Subdomain Settings . . ., and is given as 



(37) 



cbar2*x* (Haxitt*test (Haxi) 
+Hazitt*test (Hazi) 
+Hradtt*test (Hrad) ), 



(38) 



where Haxitt denotes the double partial derivative of Haxi 
with respect to time, etc. The remaining 2nd and 3rd slots of 
the dweak-tabbed panel should be zero-filled. 



18 Note that, when typing (or 'cutting-and-pasting') this and the following 
(d)weak-form expressions into their slots, all spaces and new lines must be 
eliminated from the whole expression within each slot -otherwise COMSOL 
will reject the expression. 



E. Boundary conditions 

Here the constraints stated in subsection II-C are ex- 
pressed in COMSOL-ready forms. The model resonator's 
boundary conditions are defined through the Physics =^> 
Boundary Settings . . . control panel. On the left of 
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this panel, select the Groups tab. Each named bound- 
ary Group here corresponds to a particular electromagnetic 
boundary condition, the most essential of which are described 
here. These different e.m. boundary conditions are specified 
by the expressions that populate the three slots within their 
respective contr-tabbed sub-panels, located on the right-hand 
side of Boundary Settings . . .; 'contr' here stands 
for 'constraint'. [The neighboring weak-tabbed and dweak- 
tabbed panels within the Boundary Settings . . . need 
not be touched (and left zero-filled).] 

1) Electric wall (for a bounded isotropic dielectric): 

Hrad*nx+Haxi*ny; (39) 

-Haxix+Hrady; (40) 

(Hazi*nx-Hrad*M*nx 
-Haxi*M*ny+Hazix*nx*x (41) 
+Haziy*ny*x) /x; 

here nx and ny are, as 'geometric variables' within COMSOL 
(in 2D), the components of the (outward) unit normal vector 
on the boundary of a subdomain. 

2) Magnetic wall (for a bounded isotropic dielectric): 

Haxi*nx-Hrad*ny; (42) 

Hazi; (43) 

(Haxi*M*nx+Hazi*ny 
-Hrad*M*ny-Haziy*nx*x (44) 
+Hazix*ny*x) / x. 

3) Radiation match (in free-space): As has already been 
discussed in subsection II-C, the constraints appropriate to 
implementing a radiation match, can be regarded (com- 
plex) linear combinations or 'mixings' of pure electric- and 
magnetic- wall constraints. The first constraint mixes the 
magnetic-wall constraint 23, i.e. 43 above, with the electric- 
wall constraint 20, i.e. 41 above: 

-i*cMW*Hazi*cbar*mf 

+ cEW* (Hazi*nx-Hrad*M*nx 

(45) 

-Haxi*M*ny+Hazix*nx*x 
+Haziy*ny*x) /x; 

note that 'i' here is the square root of minus one. And the 
second constraint mixes the electric-wall constraint 18, i.e. 40 
above, with the magnetic-wall constraint 22, i.e. 42 above: 

-i*cEW* (-Haxix+Hrady) 

(46) 

+cMW*cbar*mf * (Haxi*nx-Hrad*ny) . 

Here, the pair of constants {cMW and cEW}, are defined 
through TABLE III. When they are set to their standard 
(default) values of { l/\/2, 1/V2 }, equations 45 and 46 
impose a radiation match on tangential field components 
at the impedance of plane e.m. waves in free-space. Here 
also, cbar= c = 2tt/c; and mf is the mode's (center) 
frequency; both need to be defined within Options => 
Constants .... -as per their corresponding rows in TA- 
BLE III; 



The final (optional) constraint mixes the electric wall con- 
straint 17, i.e. 42 with the magnetic wall constraint 21, i.e. 44: 

tngM*cbar*mf * (Hrad*nx+Haxi*ny ) 
-tngE* (Haxi*M*nx+Hazi*ny 
-Hrad*M*ny-Haziy*nx*x 
+Hazix*ny*x) /x. 

Here; the constants tngM and tngE, are also defined through 
TABLE III and thereupon Options => Constants 
Setting {tngM, tngE} = { 1,0}({ 1,0}) constrains the magnetic 
(electric) field to be wholly tangential on the impedance- 
matching plane as is characteristic of electromagnetic trav- 
eling waves. The default setting for this third constraint was 
(arbitrarily) {tngM, tngE} = {1,0}). [It is remarked here 
that the author sought to implement the radiation-matching 
constraints more directly and elegantly with time derivatives, 
i.e., replacing 2*pi*mf *Hazi by Hazit, and similarly for 
Hrad and Haxi. But COMSOL did not generate the intended 
frequency factor when interpreting them. He thus resorted to 
entering the expressions as stated in equations 45 through 47, 
requiring mf to be set, by hand, for each mode.] 

F. Geometry 

Each resonator's geometry needs to be either constructed 
within or imported into COMSOL. COMSOL's manual pro- 
vides instructions on how to implement both. Though simple 
geometries (e.g. a cylinder of solid dielectric material inside 
a cylindrical metal can) can be quickly constructed by hand 
within COMSOL, the author found it advantageous to define 
the sets of quadrilateral subdomains into which many axisym- 
metric dielectric resonators can be readily decomposed using 
MATLAB scripts, where the script was run to generate the 
resonator's medial cross-section. The key lines in these scripts 
were those of the form 

ql = poly2vert ( [ [xl, yl] ; 
[x2,y2] ; [x3,y3] ; [x4,y4] ] ) ; 

this particular line defines a quadrilateral, named ql, 
whose vertices have the x-y coordinates: [xl,yl] , ... 
[x4,y4] . These quadrilaterals could then be imported into 
COMSOL by entering a list comprising their names into 
File => Import =^> Geometry Objects. Each complete 
MATLAB script (also available from the author upon request) 
included, where known/relevant, the cryogenic shrinkages of 
the resonator's constituent materials. 

G. Meshing 

If constructed out of quadrilaterals, the resonator's geometry 
can be meshed either into sub-quadrilaterals using Mesh => 
Map Mesh. Else, the geometry can always be meshed into 
(psuedo-random) triangles using Mesh => Initial Mesh 
with (recommended) mesh refinement over selected areas 
(covering the bright spots of WG modes). Note that, for the 
geometry to be Map Mesh-able, the vertices of its inter- 
nal quadrilaterals should generally all meet at 'cross-roads', 
where, at each, four vertices belonging to four separate quadri- 
laterals all meet at a point, as opposed to 'T-junctions', where, 
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at each, two vertices belonging to two separate quadrilaterals 
meet at a point along the boundary edge of a third quadrilat- 
eral. The reader is advised to consult the COMSOL manual 
(chapter 'Meshing', section 'Generating Meshes', subsection 
'Creating Mapped Meshes in 2D') for a fuller (though a 
still not wholly satisfactory) explanation of this rather quirky 
requirement. The meshing density can be controlled by, for 
selected edges, activating and entering an appropriate integer 
in the 'Contrained edge element distribution' 
field within of the Boundary tab of the Map Mesh control 
panel. With regard to both Map Mesh-ability and control over 
meshing density, it can be advantageous (or plain necessary) 
to divide odd-shaped subdomains (e.g. an 'L' -shaped region 
covered throughout by a single spatially uniform dielectric ma- 
terial) into several, more simply shaped adjoining subdomains 
(e.g., in the case of dividing up the L-shaped subdomain, two 
or even three rectangular subdomains). Note that the geometry 
will not mesh if the allocations of elements along too many 
edges are (inconsistently) specified; in other words, the edge 
element distribution must not be over constrained. 

H. Assignments 

1) Interiors of subdomains: The (either hand-made or 
imported) quadrilaterals composing the resonator's cross- 
sectional geometry are assigned to one of the defined 
dielectric Groups via the Subdomain Settings . . . 
=> Subdomains tab. Activating (i.e. ticking) the 
Select by group option here aids the verification 
of assignments. 

2) Edges of subdomains: The external edges of quadri- 
laterals should be set be obey one of the three above- 
described Group boundary conditions via the Physics 
=> Subdomain Settings ... =4> Subdomains tab. 
COMSOL appears to be smart enough to recognize those 
edges that are external on its own accord and automatically 
'ghosts' (grays) internal edges; the latter should not be as- 
signed to any boundary Group condition. Thus, with the 
Select by group feature activated, all the appropriate 
edges can be assigned to the appropriate (usually electric-wall) 
boundary condition in a single selection of the Groups drop- 
down menu. 

/. Solution 

In Solve => Solve Parameters: 

(a) set the selected Solver: to 'Eigenvalue:'; 
and, with the General tab selected, 

(b) set the Desired number of eigenvalues: to 
'10' -or whatever ones desires; 

(c) set Search for eigenvalue around : to '0' -or 
whatever; 

(d) set Linear system solver: to 
'Direct (SPOOLES) (this is at least the author's 
starting recommendation); 

(e) set Matrix symmetry: to Symmetric'. 
Having implemented all of the above, one should now be 

able to Solve =4> Solve Problem. 



J. Postprocessing 

COMSOL Multiphysics' standard documentation explains 
how to configure and use of the Postprocessing => 
Plot Parameters control panel. Only a few specific 
pointers are supplied here: 

[1] The center frequencies of solved resonances can be 
viewed through the Solution to use ^ Eigenvalue: 
drop-down menu in the General-tabbed sub-panel of the 
Plot Parameters control panel. 

[2] To display the morphology and features of the solved 
eigenmodes, the Expression: slot within the Surface- 
tabbed sub-panel of the same is filled with either (i) some 
function of the solved field variables {Hrad, Hazi, Haxi } or 
(ii) one of those expressions (e.g. ElecMagSqrd) pre-defined 
in Options => Expressions Scalar Expressions 
as COMSOL's equivalent of TABLE IV, or (iii) some func- 
tion/combination thereof. For example, 

loglO (AbsMagEnDens + 10~ (-10) ) (48) 

can be inserted to view the magnetic energy density on a 
logarithmic scale. To view (as a diagnostic) the divergence 
of the magnetic field strength, which should be zero, one 
inserts DivH = (Hrad-Hazi*M+ (Haxiy+Hradx) *x) /x 
instead. 

[3] Determinations of an electromagnetic mode's vol- 
ume, filling factor(s), and length (A), as per equa- 
tions 24, 25, and 27, respectively, all make use of the 
Postprocessing => Domain Integration panel. For 
example, the numerator J J J h _ s e|E| 2 dV on the right-hand- 
side of 24 can be evaluated by inserting ElecEnDens 
= Erad*Drad+Eazi*Dazi+Eaxi*Daxi into this panel's 
Expression: slot, with those entries selected in the 
Subdomain selection list on the left-hand side of the 
same panel covering all significant parts of the mode's bright 
spot. 

[4] With regards to determining filling factors, the numer- 
ator on the right-hand side of equation 25 can be evaluated 
by selecting only those subdomains filled with the relevant 
dielectric (as opposed to free-space). 

[5] To determine the resistive-wall-loss integral J |n x 
H| 2 dS (forming the denominator of equation 27), one uses 
the Postprocessing => Boundary Integration 
panel with 2*pi*x* (abs (Hazi) "2 +abs ( (Haxi*nx 
-Hrad*ny ) ) " 2 ) inserted into the Expression: slot 
within the Expression to integrate box therein, and 
where entries selected with Boundary selection cor- 
respond to the resonator's enclosing (metallic and lossy) 
surfaces. 

[6] The maximum/minimum of a field variable, as required 
to evaluate the denominator on the right-hand side of 
24, can be determined through the Postprocessing 
=> Plot parameters => Max/Min, wherein the 
Expression: slot is filled with the field variable's 
expression (viz. ElecEnDens for evaluating said 
denominator). 
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TABLE III 

COMSOL Constants -included in weak-form expressions 



1> dillC 




Expression 
(= value) 




uesciipuon [ uiiii j 


C 




299792458 




~a — FT — u — 7 w — Ti — 

speed or light (exact) [m/s] 


cbar 




2*pi/ c 




frequency constant [s/mj 






( = 2.095845e-8) 




cbar 2 




4*pi " 2/ c" 


2 


2 2 

frequency constant [s /m J 






( = 4.392566e-16) 




alpha 




1.0 




penalty-term coefficient 


M 




9 




azimuthal mode order 


eO 




1.0 




relative permittivity 










of free-space 


el 




n_AlGaAs" 


2 


relative permittivity of 






( = 11.2896) 




isotropic_dielectric_ 1 


e2 




1.0 




same but of 










isotropic dielec._2 


e3 




1.0 




etc. 


eperpO 




1.0 




relative permittivity 










oi tree-space 










in directions transverse 










to cylindrical axis 


eparaO 




1.0 




same but in direction 










parallel to cylindrical axis 


eperpl 




9.2725 




relative permittivity of 










uniaxial_dielectric_ 1 










in directions transverse 










to cylindrical axis 


eparal 




11.3486 




same but in direction 










parallel to cylindrical axis 


epara2 




1.0 




relative permittivity of 










uniaxial_dielectric_2 










transverse to cyl. axis 


eperp2 




1.0 




same but parallel to 










cylindrical axis 


eperp3 




1.0 




etc. 


epara3 




1.0 




etc. 


e_2 93K_alumina 


9.8 




relative permittivity of 










alumina at room temperature 


epe_4K_sap_ 


_UWA 


9.2725 




UWA values for cryogenic 










HcJVLbX sappnire 


epa_4K_sap_ 


_uwa 


1 1.3486 






epe_2 93K_sap 


9.407 




nominal room temperature 










values tor same 


epa_2 93K_sap 


11.62 






epe_4K_sap_ 


_NPL 


9.2848 




Values fitted to 










JNrL Cs-rountain 










HbJVlhX resonator 


epa_4K_sap_ 


_NPL 


11.3660 






n_silica 




1.4457 




refractive index of 










thermally grown 










silica (Fig B.2, p. 172 of 










ref. [28]) 


n_AlGaAs 




3.36 




average refractive index of 










GaAs and AlGaAs layers 










(p. 1/2 or ret. |_14J) 


mf 




2.374616el4 




1 1 KILL 1 ) 11 CUUCllLV 


ttgH 




1 




toggle 


ttgE 









toggle 


mix_ang 




45 




electric-magnetic mixing 










angle (in degrees) 


cMW 




sin (mix_ang 








* pi /18 


0) 


magnetic-walled-ness 






(= 0.707107) 






cEW 




cos (mix_ang 








* pi /18 


0) 


electric-walled-ness 






(= 0.707107) 






tngM 




1 






tngE 












TABLE IV 

COMSOL Scalar Expressions -for postprocessing 



Name 


Expression 




(Description) 


DivH 


(Hrad-Hazi*M+ (Haxiy+Hradx) *x) /x 




(divergence of magnetic field -should be zero!) 


MagEnDens 


Hrad*Hrad+Hazi *Hazi + Haxi *Haxi 




(magnetic energy density) 


Drad 


(Haxi*M-Haziy*x) /x 




(radial component of electric displacement) 


Dazi 


-Haxix+Hrady 




(azimuthal component of electric displacement) 


Daxi 


(Hazi-Hrad*M+Hazix*x) /x 




(axial component of electric displacement) 


Erad 


Drad/eperp 




(radial component of electric field strength) 


Eazi 


Dazi/eperp 




(azimuthal component of electric field strength) 


Eaxi 


Daxi/epara 




(axial component of electric field strength) 


ElecMagSqrd 


Erad*Erad+Eazi *Eazi+Eaxi*Eaxi 




(electric field strength magnitude squared) 


ElecEnDens 


Erad*Drad+Eazi*Dazi+Eaxi*Daxi 




(electric energy density) 


AbsMagEnDens 


abs (Hrad) "2+abs (Hazi) "2 




+abs (Haxi) "2 




(absolute magnitude energy density) 


MagNrmlHSqrci 


2*pi*x*abs (Haxi*ny 




+Hrad*nx) "2 




(magnitude normal mag. field strength squared) 


MagTngHSqrd 


2*pi*x* (l*abs (Hazi) "2 




+l*abs (Haxi*nx-Hrad*ny) " 2) 




(magnitude tangential magnetic field squared) 


AbsElecSqrd 


abs (Erad) "2+abs (Eazi) "2 




+abs (Eaxi ) " 2 




(absolute electric field squared) 



